#ifndef _FLOWSIM_HEADER_
#define _FLOWSIM_HEADER_

#include <memory>
#include <vector>
using namespace std;

#include "Ogre.h"
#include "FlowSimHeaders.h"
#include "VectorField2.h"
#include "VectorField2Render.h"
#include "ScalarField2.h"
#include "DyeField3.h"
#include "TxtDebugOut.h"
#include "DyeFieldRender.h"

namespace FlowSim
{
  /**
   *This class is responsible for the simulation of an incompressible flow, governed by the the Navier-Stokes
   *equations for incompressible flows. The simulation is stable, per the method seen in Jos Stam's 1999 paper: 
   *Stable Fluids. We define the simulation as:
   *S(w) = P o F o D o A(w), a composition of a set of transformation on a vector field u.
   *Where,
   *A(w) : Advection.
   *D(w) : Diffusion.
   *F(w) : Force
   *P(w) : Projection. Where P projects a vector field w onto it's divergence-free component, utilizing the concept
   *       of the Helmholtz-Hodge decomposition, such that, w = u + Del p. Where u is divergence free, and Del p is
   *       the graident of a pressure field p. Thus, u = w - Nabla p. Where Del is the Nabla operator.
   **/
class MyFlowSim
{
 
public:
  MyFlowSim();
  ~MyFlowSim();
  void init(Ogre::SceneManager* sMgr);
  void update(); //update the simulation
  void setInput(VectorField2* vf); //input: initial vector field.
  void testInterpolation(int x, int y, int wWidth, int wHeight);
  Ogre::Real getDT(){return _dt;}
  void setDT(Ogre::Real dt){_dt = dt;}
   //testing methods
  void testAdvectV();
  //reset
  void reset();
  void addSource(int x, int y, int wWidth, int wHeight);
   void addForce(int x, int y,int width,int height);
  void moveLeft();
  void moveRight();
  void moveUp();
  void moveDown();
  void changeSpeed(int speed);
  void changeDifRate(float dr);
  void toggleVFRender();
protected:  
  //The fields.
  auto_ptr<VectorField2> _vf; //the velocity field
  auto_ptr<VectorField2> _tvf; //temp velocity field
  auto_ptr<ScalarField2> _psf; //pressure scalar field.
  auto_ptr<ScalarField2> _tpsf; //temp pressure scalar field.
  auto_ptr<ScalarField2> _divwsf; //divergence of w scalar field.
  auto_ptr<DyeField3> _dyef; //the dye field.
  auto_ptr<DyeField3> _tdyef; //the temp dyefield.
  //dye source vector
  vector<int> _vsx;
  vector<int> _vsy;
  //renderers
  VectorField2Render* _vfRender;
  DyeFieldRender* _dyeRender;
  //simulation variables
  Ogre::Real _dt; //delta time for the flow simulation.
  int _dim; //simulation dimension
  Ogre::Real _vis; //simulation viscosity variable
  Ogre::Real _sgr; //grid resolution, i.e the dimension between each grid cells.
  int _wallSize; //wall size
  Ogre::Real _dRate; //rate of diffusion
  Ogre::Real _windSpeed;
  Ogre::Vector2 _windForce;
  //misc
  TxtDebugOut* _dout;
  bool _brVF;
  //graphics
  Ogre::SceneManager* _scnMgr;
  int _playerX,_playerY;
  int _playerDim;

  //Methods:
  void advectVel(auto_ptr<VectorField2> &vf,auto_ptr<VectorField2> &temp); //advection of velocity
  void boundaryVel(auto_ptr<VectorField2> &vf); //handle velocity boundary conditions
  void boundaryPre(auto_ptr<ScalarField2> &sf); //handle pressure boundary conditions
  void boundaryDye(auto_ptr<DyeField3> &df); //handle dye boundary conditions: NOT USED.
  void project(); //project to get divergence free velocity field
  //helper methods
  Ogre::Vector2 lerp(Ogre::Vector2 z1, Ogre::Vector2 z2,
		     Ogre::Real x1, Ogre::Real x2,Ogre::Real x); //linear interpolation, Vector 2
  void lerp(Ogre::Real y1r, Ogre::Real y1g, Ogre::Real y1b,
	    Ogre::Real y2r, Ogre::Real y2g, Ogre::Real y2b,
	    Ogre::Real x1, Ogre::Real x2, Ogre::Real x,
	    Ogre::Real &rr, Ogre::Real &rg, Ogre::Real &rb); //linear interpolation, uint8
  void bilerpDye(const Ogre::Vector2 p, Ogre::Real &r, Ogre::Real &g, Ogre::Real &b); //bilerp for dye.
  Ogre::Vector2 bilerp(auto_ptr<VectorField2> &vf, const Ogre::Vector2 p); //bilerp for velocity, vector2
  void swap(auto_ptr<VectorField2> &v1, auto_ptr<VectorField2> &v2); //swap vector fields
  void swap(auto_ptr<ScalarField2> &v1, auto_ptr<ScalarField2> &v2); //swap scalar fields
  void swap(auto_ptr<DyeField3> &df1, auto_ptr<DyeField3> &df2); //swap dye fields
  void jacobi(auto_ptr<ScalarField2> &x, auto_ptr<ScalarField2> &tempx, 
	      auto_ptr<ScalarField2> &b,Ogre::Real alpha,Ogre::Real rBeta,int numIt); //the jacobi iteration method with scalar fields x,b.
  void jacobi(auto_ptr<VectorField2> &x, auto_ptr<VectorField2> &tempx,
	      auto_ptr<VectorField2> &b, Ogre::Real alpha, Ogre::Real rBeta, int numIt);
  void divergence(auto_ptr<VectorField2> &vf, auto_ptr<ScalarField2> &sf); //compute the divergence of a vector field.
  void computeDivFreeField(auto_ptr<VectorField2> &w, auto_ptr<ScalarField2> &p,
			   Ogre::Real halfrdx); //compute divergence free vector field.
  void sourceDye(auto_ptr<DyeField3> &df); //add source to dye field
  void addConstVel(auto_ptr<VectorField2> &vf);
 
  void advectDye(auto_ptr<VectorField2> &vf, auto_ptr<DyeField3> &df,
		 auto_ptr<DyeField3> &tdf);
  void testAdvectVel(auto_ptr<VectorField2> &vf,auto_ptr<VectorField2> &temp);
  void resetSink(auto_ptr<VectorField2> vf,
		 auto_ptr<VectorField2> tvf,
		 auto_ptr<ScalarField2> psf,
		 auto_ptr<ScalarField2> tpsf,
		 auto_ptr<ScalarField2> divwsf,
		 auto_ptr<DyeField3> dyef,
		 auto_ptr<DyeField3> tdyef); //reset sink
  void resetSink(auto_ptr<VectorField2> vf,
		 auto_ptr<VectorField2> tvf,
		 auto_ptr<ScalarField2> psf,
		 auto_ptr<ScalarField2> tpsf,
		 auto_ptr<ScalarField2> divwsf);
	
 
 //void setRightBoundVel(); 
};
}
#endif
